A new effective interaction for the trapped fermi gas: the BEC-BCS crossover 
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We extend a recently introduced separable interaction for the unitary trapped Fermi gas to all 
values of the scattering length. We derive closed expressions for the interaction matrix elements and 
the two-particle eigenvectors and analytically demonstrate the convergence of this interaction to the 
zero-range two-body pseudopotential for s-wave scattering. We apply this effective interaction to the 
three- and four-particle systems along the BEC-BCS crossover, and find that their low-lying energies 
exhibit convergence in the regularization parameter that is much faster than for the conventional 
renormalized contact interaction. We find similar convergence properties of the three-particle free 
energy at unitarity. 
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I. INTRODUCTION 

Atomic Fermi gas systems have generated much in- 
terest in the last several years [3, 0] • While these sys- 
tems can be well-described by a relatively simple model 
Hamiltonian, at low temperatures they exhibit a rich phe- 
nomenology whose theoretical description has proven to 
be challenging. At low momenta, or when the range of 
the interaction is sufficiently short, the dominant scat- 
tering process occurs in the s-wave channel, so that a 
single parameter, the s-wave scattering length a, suffices 
to characterize the inter-particle interaction. Depend- 
ing on the value of fcp a (where fcf is the Fermi momen- 
tum), widely different behavior is observed. This ranges 
from the formation of tightly bound dimers in the Bosc- 
Einstcin condensate (BEC) regime when kpa is small 
and positive, to Cooper pairing in the Bardeen-Coopcr- 
Schricffcr (BCS) regime when kpa is small and negative. 
These behaviors are connected by a strongly interacting, 
nonperturbative regime, known as the unitary regime, 
where \a\ is much larger than any other length scale in 
the system. Remarkably, each of these regimes is acces- 
sible experimentally, and all exhibit superfluid behavior 
below a certain a-dependent critical temperature. 

While accurate mean-field theories exist for the BEC 
and BCS regimes, there is no simple approximation 
that can accurately describe the transition between these 
regimes through the unitary limit. As a result much in- 
terest has been taken in applying numerical methods, 
such as quantum Monte Carlo and numerical diagonal- 
ization. For systems varying from ~ 10 to hundreds of 
particles, quantum Monte Carlo calculations have been 
carried out to calculateground state P4l3j and ther- 
modynamic properties |M| - tl8| . For smaller numbers of 
trapped atoms, energy spectra have been calculated using 
a basis set expansion method with correlated Gaussians 
in coordinate space (for up to six atoms) [1, |H| and a 
stochastic variational approach (for four atoms) 19[ . 

The three-body problem has been solved analytically 
in the unitary regime [20l l2~l| and numerically to high 
accuracy along the BEC-BCS crossover (22), [23[- In fact 
the few-body trapped cold atom problem has become 



more interesting recently, since it was pointed out that 
the virial expansion for the partition function in a har- 
monic trap works well at unexpectedl y lo w temperatures 
into the quantum degenerate regime [231 ] . Moreover, the 
scaling relations which exist for the dilute gas can allow 
larger systems to be addressed from the study of smaller 
ones HIS!- 

Here we discuss the cold atom problem in the con- 
text of the configuration- interaction (CI) approach to in- 
teracting many-particle fermionic systems. In this ap- 
proach, widely used in atomic, molecular and nuclear 
physics, a single-particle basis in the laboratory frame 
is used to construct a many-particle basis of Slater de- 
terminants for a fixed number of fermions. The Hamilto- 
nian matrix is then calculated in this many-particle basis 
and diagonalized. For a harmonically trapped system, a 
natural choice for the single-particle basis is that of the 
three-dimensional harmonic oscillator. 

The short-range interaction of the cold atom problem 
is often approximated by a zero-range interaction. How- 
ever, the obvious form of such a potential, a contact inter- 
action Vo8(r), is ill-defined in three dimensions and must 
be regularized. This is usually accomplished by intro- 
ducing a momentum or energy cutoff in relative motion 
and renormalizing the strength Vq of the interaction so it 
reproduces the physical scattering length for the uniform 
gas, or the lowest bound state energy in relative motion 
for the trapped system. The many-particle energies are 
then calculated as a function of the regularization cutoff 
parameter. 

A contact interaction is non-vanishing only for relative 
angular momenta / = of the two particles. For a basis 
of harmonic oscillator wave functions, a natural regular- 
ization parameter is the number of oscillator s waves in 
relative motion. However, the convergence of the many- 
particle energies as a function of this regularization pa- 
rameter for the renormalized contact interaction is slow 
(as a low negative power) (2(| [27} . 

In Ref . [27] a new effective interaction was introduced in 
the unitary limit for which the low-lying energies of the 
three- and four-particle systems converge substantially 
faster than for the conventional renormalized contact in- 
teraction as a function of the regularization parameter. 
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While this effective interaction is no longer a contact in- 
teraction, it reproduces the same many-particle energies 
in the limit when the regularization parameter is sent to 
infinity. The faster convergence of the many-particle en- 
ergies enables their calculation to higher accuracy when 
using this new effective interaction. Another approach 
to improve the convergence relative to the renormalizcd 
contact interaction was investigated within the frame- 
work of effective field theory by including perturbatively 
next-to- leading order and next-to-next-to-leading order 
interactions [28|. 

Here we generalize the construction of such an effec- 
tive interaction away from unitarity. We study the low- 
energy spectra of the three- and four-particle systems and 
find convergence properties versus the regularization pa- 
rameter along the complete BEC-BCS crossover that are 
similar to the unitary regime. We also use this effective 
interaction to study the thermodynamics of the three- 
particle system at unitarity. In particular, we demon- 
strate the exponential convergence of the free energy at 
finite temperature as a function of the regularization pa- 
rameter. The converged free energy is compared with the 
exact free energy constructed from the known analytical 
spectrum of the three-particle system [2(| Hl| . 

This paper is organized as follows. In Sec. II, we briefly 
review the two-particle cold atom problem in a harmonic 
trap. In Sec. Ill, we derive closed expressions for the ma- 
trix elements of the new effective interaction for any scat- 
tering length and demonstrate convergence (as a func- 
tion of the regularization parameter) of its two-particle 
eigentates to the exact eigenstates in the unitary limit. 
In Sec. IV, we use this interaction to study the conver- 
gence (vs. the regularization parameter) of the spectra 
of the three- and four-body systems along the BEC-BCS 
crossover and the three-particle free energy at unitarity. 
The latter is compared with exact results. Finally, in 
Sec. V we present our conclusion. 



where g — 2Trfi 2 a/fi relates the scattering length a to 
the strength of the interaction Q ■ Equivalently, one can 
impose the Bethe-Peierls contact conditions on the wave- 
functions; see, e.g., in Ref. HtJ Several interesting con- 
sequences of the form ([2]) of the interaction potential for 
the many-body system have been explored in Ref. l29l 

For a two-particle system, the Hamiltonian is separable 
in the center-of-mass and relative coordinates R = (ri + 
r2)/2 and r = — r\. so that H = Hqm + H Te \, with 



H Ie 



2/z 



-AiwV + V(r) = H + V(r) 



where Hq is the harmonic oscillator Hamiltonian in rela- 
tive motion with reduced mass fi = to/ 2. 

The eigenstates of the non- interacting two-particle sys- 
tem in a harmonic trap may be labeled as \AfCA4nlm), 
where J\f, C, M. are the radial, angular momentum, and 
magnetic quantum numbers for the center of mass mo- 
tion, and n,l,m are the corresponding quantum num- 
bers for the relative motion. The associated energies are 
E = (27V + £ + 3/2 + 2n + / + 3/2)fra;. A pure s-wave in- 
teraction leaves the I ^ states and energies unchanged 
while mixing the I = states into eigenstates we denote 
by \J\f£Mu®), with E = 2A/'+£ + 3/2 + £ i . The energies 
Si in relative motion for scattering length a are give n by 
the solutions to the transcendental equation j3oll3l1] 



r(-g/(2fta;) + 3/4) _ 1_ 

T(-e/(2fkJ) + 1/4) ~ V2a/a 



(3) 



The exact wave functions in relative motion arc also 



known. They arc given by [31 



^> = £4 i) |n00>, 



71=0 



II. TWO-PARTICLE PROBLEM 

The trapped two species cold atom system is modeled 
by the Hamiltonian 



N fc2 » , 



(i) 



i=l 



i=l 



where N — N\ + N2 , -/V 1 and N2 are the number of atoms 
for each species (spin- up and spin-down), lo is the trap 
frequency, and V{r) is a short-range interaction. A nat- 
ural length scale for this Hamiltonian is the oscillator 
length a osc = ^Jh/muj. 

The interaction V(r) is modeled by a zero-range, pure 
s-wave interaction, which due to Pauli exclusion acts only 
between particles of differing species. This can be ex- 
pressed by a regularized (5-function, 



V(r) = gS 3 (r)(d/dr)r , 



(2) 



where 



<(0) 



Here a n = (2n + 3/2)hui are the non-interacting rela- 
tive energies, Ai is a normalization factor, and <p n (r) = 
(Pnoo(r) is a harmonic oscillator wave function, with 



^„(o) = ^- 3 / 4 



(2n + l)!! 



(2n) 



1/2 



In the unitary limit, where e% = (2i + \/2)hw, the 
normalization factor Aj has the simple form 



A: 



_ T -3/2^ (2"+l)!! 



1 



(2n)!! [2{n-i) + l] 2 



TT- 1 / 2 (2i)\\ 
2 {2i-\)\\ 



3 



III. EFFECTIVE INTERACTION 

Wc first discuss the conventional contact interaction. 
It can be regularized by introducing an energy cutoff in 
relative motion. This defines a sequence of interactions 
V c iq \ q = 0,1,2,..., with 

,0^m,0 , (4) 

for n, n! = 0,1,..., q, and where tp n (0) = [(2n + 
l)!!/(2?i)!!] 1/2 . All other matrix elements of V c {q) are 
taken to be zero. The coefficient \g is determined for 
each q so as to reproduce the lowest relative-motion en- 
ergy £o. We refer to this interaction as a renormalized 
contact interaction. 

The contact interaction in Eq. (|4| is separable in the 
relative oscillator basis. A general separable interaction 
for / = has the form 

(nlm\V$\n'lm) = f*f n >8i, 5 m ,o (n,n'<g). (5) 

An effective interaction for the trapped cold atom sys- 
tem can be determined by choosing the coefficients /„ to 
reproduce the lowest q+1 relative energies eq,... , £q of 
the two-particle system. This idea was used in Ref. [27| to 
construct a new effective interaction in the unitary limit 
of infinite scattering length. Here we consider the more 
general problem for any value of the scattering length a. 



A. General problem and solution 

The matrix elements of the relative-coordinate Hamil- 
tonian take the form 

(n00|#<$ |n'00) - 6 n , n ,a n - /*/„, , (6) 

where a n are the non-interacting eigenvalues of Hq. To 
determine the coefficients /„ we derive the following re- 
sult. 

Theorem. Let eq, e\, . . . , e q and ao, ct\, . . . a q be real 
numbers such that Eq < ... < e q , ao < ... < a q and 
Ei 7^ OLj for all i,j, and let fa, fx, ■ ■ ■ , f q be complex. Then 
the (g + l)-dimensional matrix H n ^ n > = a n S n ,n' — f„fn' 
has the Ei as its eigenvalues if and only if 



Hk( a n-£k) 



V llfc#n( a « -«fc)' 

in which case the ith eigenvector feW has components 



(7) 



fn 



a n - £i 



C, 



ILK - £ ^ 



(8) 



A solution satisfying (J7J exists if and only if £o < ao < 
E\< ... <a q . 



Proof. An eigenvector &W with eigenvalue Ei must sat- 
isfy 



n=0 



which implies 



Otn - Ei 



(9) 



(10) 



where Cj = ^TJ m f m bm is a normalization constant. In- 
serting (fTU|) into ©, we obtain 



<i 

E 

n=0 



| ^ 



where the matrix M is defined by 
M in = 



(11) 



Thus the f„ are determined from the eigenvalues by the 
matrix equation Mv = 1, where 1 = (1,1,...,1) and 

« = (l/olM/i| 3 ,...,l/ 9 IT- 

The Matrix M is of a well-known form, called a Cauchy 
matrix, whose inverse is given by [32j 



(a„ - ei) il^o^ - ^) il^k - 

We therefore have 
|/ n | 2 = (M" 1 !^ 



E 



n fe ("™ ~E k ){a k -Ei) 



{a n - Ei) Uk^i( £ k ~ Uk^n( a n - a k) 

n fe (a„-e fe ) ^ I\k^J £ i - a k) 



ilynK " a k) i IL^(£i ~ £ k) 



(12) 



The sum in the last expression can be evaluated using 
the identity [H] 



0, < r < N - 1 

1, r = N-l 



E 



(13) 

where xi,...,xjv are all distinct. This is done by ob- 
serving that rife^n( £ » — a k) — £ i + P( £ i), where P is a 
polynomial of degree smaller than q. There are q + 1 of 
the £j, which we take to be Xi, X2, Xg+i in Eqs. (|13l) . 
According to the first case in Eqs. (fl3| . the contribution 
to the last sum in Eq. (|12|) from P(£i) vanishes, while ac- 
cording to the second case of (fTB")) . the contribution of e? 
to the sum is exactly 1. Thus, the components /„ satisfy 
(J?]) . On the other hand, it is easy to verify that if the /„ 
satisfy {7J) then given by (JTOJ) is an eigenvector with 
eigenvalue Ei. 
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If £o < cko < £1 < • ■ ■ < ot q then it can easily be 
seen that the argument to the square root in Eq. ([7]) is 
positive, and therefore a solution for the /„ exists. Con- 
versely, we can show that if the argument to the square 
root is positive for all n, the eigenvalues must have this 
ordering. For each n = 0, . . . , q let k n be the number of 
eigenvalues £j for which £, < a„. Then 



sign 



In order for the sign to be positive for all n, we therefore 
require that n + 1 — k n is always even. But since < 
ko < . . . < k q < q + 1, the only possibility is k n = n + 1, 
i.e., £q < «o < £i < • • • < OL q . 



The eigenvectors are given by Eq. (|10[) . To determine 

the normalization constant C{, we require X)n=o l^n^| 2 = 
1 and use Eq. (JT]) for |/„| to find 



g 



n 



(a„ - £fc) 



1 



This sum can also evaluated using the identity (|13l) . To 
do this, we define a q +i = Si and rewrite 



c; 



n 



k^i,0^k^q 



-0 llk^nfi^k^q+l^ 

We now add and subtract the n = q + 1 term to obtain 

Uk^.0^q( a n ~ Efc) 



n=0 Ylk^n,0^k^q+l( a n a k) 

Taking xi = ao,X2 = ai,...,ccjv = otq+i, the first case of 
(|13p may be applied to find that the sum vanishes and 
thus 



c; 



rifeK - £ i) 



This completes the proof. 

Since the eigenvalues £j, a„ for the cold atom problem 
are indeed ordered as the theorem requires, we choose 
fn = \fn\ to be positive real numbers according to 
Eq. ©. 

The complete set of the two-particle eigenvectors of the 
effective interaction v£§' is therefore given by (i) the non- 
interacting states \N CAinlm) for I > or (n > q and I = 
0) with eigenvalues (27V + C + 3/2 + 2n + I + 3/2) huj, 
and (ii) the I = interacting states 



\Af£Mb {i 



Uk( a k ~ gj) 
Yl k7 ti( £ k - Ei) 



1/2 



n=0 



In 



-\NCMn00) (14) 



with eigenvalues (27V + £ + 3/2 + Si) huj and where the 
£j are given by the solutions of Eq. ([3]). 

We also derive closed expressions for the trace and 
norm of the (q + l)-dimensional matrix IQj defined by 

(V}s)nn' = (nOOlV^V'00). Its trace has a simple form 
which can be obtained via Eqs. (J7J and (fT3"]) , Taking 
Xi = ctQ,X2 = a±,...,XN = Q!g and applying the second 
and third cases of Eq. (TT3"|) , we find 



n=0 

This result implies that the Frobenius norm of vjg is 



(<?)\ |2 



ynn' | I 
n.n' \n— 



)](a n - £„) 



B. Properties in the unitary limit 

In this section we consider the properties of the separa- 
ble effective interaction in the unitary limit. The simple 
form of the relative-motion unitary eigenvalues £, allows 
us to verify Eq. (9) in Ref. for the interaction pa- 
rameters /„ at unitarity and prove analytically that 
converges to the pseudopotcntial ([2]) in the limit q — > oo. 



1. Interaction parameters f n 

In the unitary limit e, = 2i + 1/2 are the harmonic os- 
cillator energies shifted down by one oscillator quantum. 
In this case the expression for /„ reduces to 



f2 _ nLo(2("-fc) + i) 

Jn Uk^(n~k) 

IlLo(2fc + l)lE(-2fc+l) 

_ (2n + l)!! (2(g-n)-l)!l 
~ (2n)!! (2(g-n))!! ' 

which is exactly Eq. (9) of Ref. H3. 

2. Convergence of two-particle eigenvectors 



(15) 



To prove the convergence of the effective interaction 
to the pseudopotential V(r) in Eq. ([2]), we study the 
convergence of the two-particle eigenstates \b^) of Hq + 
v}fj) to the corresponding eigenstates \u^) of Hq + V(r) 
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as q — >• oo. For any scattering length, the square of the 
n-th component of \b^) is given by 



states in the harmonic oscillator shells N = 0, 



1 



]lfc^( £ fe " £ i) Uk^n( a n ~ <*fc) (a„ - e l ) 2 ' 
In the unitary case 

(hil)) 2 = II fc (2(fc-0 + l)II fc (2(n-fc) + l) 1 
l " j n fc#i 2(A-i) n M „2(n-fc) (a n -^) 2 
1 (2i-l)!!(2( 9 -i) + l)!! 



(a„- £i ) 2 (2*)!! (2(g -»))!! 
(2n+l)!!(2(g-n)-l)H 



(2n)H (2(g-n))H ' 
The asymptotic behavior for large q is 



# } (<z) 



(a n - Si) 



(2i - 1)!! 
(2i)!! 



1/2 



(2n+l)!!l 1/2 



(2n) 



9 - i + 3/4 



g-n+1/2 



i/4 



We compare to the exact wave function \u^) by comput- 
ing the relative difference for each component, 



1 



1/4 



1/4 



q - n + 1/2 

n- i + 1/4 
~ 4(g - n + 1/2) ' 

This diminishes like l/q for large g, proving convergence 
to the exact two-body eigenstates in relative motion. Im- 
portantly, this is faster than the 1 iJ q convergence of the 
renormalized contact interaction 27 1. 



IV. APPLICATIONS 

To demonstrate the advantage of the new effective in- 
teraction as compared with the renormalized contact in- 
teraction, we study the convergence versus q of low-lying 
energies of the three and four-particles systems at var- 
ious values of the scattering length. We also study the 
convergence versus q of the three-particle free energy at 
unitarity (for which an exact solution exists). The latter 
is important for thcrmodynamical studies of the trapped 
gas. The three-particle system we study consists of two 
spin- up particles and one spin-down particle (tti)j while 
the four-particle system is unpolarizcd (TtlD- 



A. Spectroscopy 

The CI method works in the laboratory frame and re- 
quires a certain truncation scheme. As in Ref.[27l we use 
as our many-particle basis the set of all Slater determi- 
nants which can be constructed from the single-particle 



For a given regularization parameter q, the two-body in- 
teraction matrix elements in the laboratory frame are 
calculated by transforming to the relative and center of 
mass coordinates via the Talmi-Moshinsky brackets [33| 
and using the interaction ([5]) in relative motion. We carry 
out direct diagonalization of the CI Hamiltonian using a 
new code for three and four particles which works in a 
two-species formalism with a basis of good total orbital 
angular momentum L and parity tt. In order to deter- 
mine the many-particle energies, we take the following 
limits. For fixed q we calculate the energies as a func- 
tion of iVmax and extrapolate to N max — > oo following 
the method of Ref. [27]. The resulting energies are then 
studied as a function of q to obtain an estimate of the 
q — > oo limit. 

In the CI calculations we used PARPACK for paral- 
lel diagonalization on a Linux cluster using between 1 
and 320 CPU cores. The largest system studied, the 
four-particle system with A max = 12 and total L = 1 
(involving ~ 4.6 million configurations), required a to- 
tal of about 2300 GB to store the matrix elements of the 
Hamiltonian and 24 hours to run on 320 cores, while more 
typical calculations, e.g., four particles at -/V max = 9 and 
L = 1 (520, 000 configurations) required about 64 GB of 
storage and four hours of computation time on 32 cores. 



1. Three particles 

We study a selected set of the low-lying energies of the 
three-particle system as a function of the inverse scatter- 
ing length. Here the main topic of interest is the con- 
vergence of the eigenvalues (i) as the size of the model 
space (i.e., the maximal number of oscillator shells iV max ) 
increases, and (ii) as the regularization parameter q in- 
creases. 

(i) Overall, convergence in N max is fast for both the 
renormalized contact interaction and the new effective 
interaction, and similar to what was already discussed in 
Ref. [!?] in the unitary limit. Convergence slows some- 
what as we follow the crossover from the BCS to the 
BEC regime and tends to require more shells as q in- 
creases. At q = 4, N max = 13, we find that the lowest 
four eigenvalues for each of L x = + , 1 + , 1~, 2 + , and 
2~ are converged in N max to about 0.0003% in the BCS 
regime (a osc /a = —0.8), 0.002% in the unitary limit and 
0.01% in the BEC regime (a osc /a) = 0.8), on average. 

(ii) As a function of q (for 1 ^ q ^ 6), we iden- 
tify two convergence patterns for low-lying three-body 
eigenvalues. The first occurs for most eigenvalues, and 
is characterized by fast, monotonic convergence for 1 

g ^ 4 or 1 ^ g ^ 5 at an exponential or faster-than- 
cxponential rate. The values at larger q then give small, 
often nonmonotonic corrections. As an example, the 
lowest L n = + state at unitarity has the exact value 
E = 4.6662 . . . hot. With the effective interaction, we 
find the energy of this state decreases monotonically to 
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FIG. 1. Convergence of the lowest 1- state for the three parti- 
cle system on either side of the crossover. The energy IS 
plotted as a function of the regularization parameter q. Solid 
circles are for the effective interaction and open circles are 
for the renormalized contact interaction, (a): The BCS side, 
a OS c/a = -0.8. (b): The BEC side, a OBC /a = +0.8. Inset: 
the differences AE {q) = E {q) - E ( -"- 1) on a log scale. For 
the effective interaction the energy decreases slightly above 
q = 4, with the q = 6 values equal to 4.7392 fku and 3.2664 fiu 
in cases (a) and (b), respectively. 




FIG. 2. (Color online) Selected low-lying energies of the three- 
particle system vs. scattering length. Symbols are the esti- 
mates for q — > oo as described in the text. Lines are linear 
interpolations. 



E = 4.6602 hw at q — 4 (slightly overshooting), and 
E = 4.6597 hui at q = 5. The convergence then switches 
direction at q = 6 and approaches the exact value with 
a low power-law behavior, so that the q 6 values pro- 
vide only a 0.14% correction to the q — 4 value. For 
such eigenvalues we report the energy at the highest q 
calculated, usually q — 6. 



The second convergence pattern we encounter is com- 
pletely monotonic for 2 ^ q ^ 6, but with a noticeably 
slower rate of convergence. For these eigenvalues, the 
convergence is smooth, and we estimate the q — > oo value 
by fitting \og(AE^), where AE^ = E^ - E^' x \ to 
a quadratic polynomial. This provides an interpolating 
function f(q) = A + Bq + Cq 2 (with C < 0) which we 
then use to numerically sum the residual terms in the 
expansion £(°°) = + A£^ +1 ) + AE^+ 2 ^ + ... m 
£(<?) + 10/(9+1) + l0/(5+ 2 ) + . . .. An example of such a 
state is the 0^" state at a osc /a = 0.0. (We label energy 
levels as L™, where n indexes the levels with total angular 
momentum L and parity ir in increasing order starting 
from zero, so that 0^" is the fourth-largest eigenvalue of 
those with L = 3 and positive parity.) In this case, the 
q = 4 value is 6.716 hw and the extrapolated value (from 
q = 6) is 6.657, indicating the state with the exact energy 
of 6.6662 ...hui,& 0.14% error. 



Comparing with the renormalized contact interaction, 
we find significantly faster convergence with the new in- 
teraction for both convergence patterns. For example, 
the 0^" state at unitarity attains the value 6.984 hui at 
q = 4 with the contact interaction, a difference of 5% 
from the correct value, compared to 0.7% at q — 4 for 
the new interaction. The difference for faster-converging 
eigenvalues is not as dramatic, but the new interaction is 
still clearly preferred. This is shown in Fig. Q] for the 1q 
state on either side of the crossover, at a osc /a = —0.8 and 
Oosc/o = 0-8- Importantly, we see that the advantage of 
the new interaction is maintained away from unitarity. 
This holds true for all the eigenvalues we studied, and 
for the four-particle system as well. 



Because convergence is only truly attained for large q, 
and it is difficult to assess the effects of non-monotonicity, 
we estimate the systematic error in our three-particle en- 
ergies along the crossover based on comparison with the 
unitary calculations, where we compare the convergence 
in q with the known exact results. From this we esti- 
mate an upper bound for the systematic error due to non- 
convergence of approximately 0.2% in our three-particle 
energies. 



Using the new effective interaction, we also present in 
Fig. rj] results for a selected set of the low- lying energies 
of the three-particle system as a function of the inverse 
scattering length. Although the three-particle system is 
now well-studied, this provides a clear verification of the 
capability of the method. Our results agree with those 
published previously [U H(| . As in the the earlier stud- 
ies, we find crossings between levels with different quan- 
tum numbers, for instance the 2^ state and the 1]~ state 
at approximately a osc /a = 0.14, and the 1q state and the 
0q state at approximately a osc /a = 0.93, which changes 
the parity and angular momentum of the ground state. 



7 



2. Four particles 



The four-particle system is less well-studied and pro- 
vides a more interesting test of the new effective interac- 
tion. As with the three-particle system, we examined a 
selected set of the low-lying energies as a function of the 
inverse scattering length. 

(i) Overall we find fast convergence in N maXl very sim- 
ilar to the three-particle case. More precisely, at q — 4, 
Nmux = 11, the level of convergence ranges from an 
average of about 0.05% at a osc /a = —0.8 to 0.2% at 
Oosc/ a = 0.8 for the lowest four states of each of the 
0~, + , 1~ and 1 + configurations. 

(ii) As a function of q, the two convergence patterns ob- 
served for the three particle case - the more common fast, 
slightly non-monotonic convergence and the less common 
slower, monotonic convergence - are still prevalent, but 
with the addition of a few eigenvalues which are non- 
monotonic for smaller q. An example is the 0^~ state at 
unitarity, which becomes non-monotonic at q = 3. For 
eigenvalues with the first two convergence patterns, we 
estimate the q — > oo values in the same manner as the 
three-particle case. For the new non-monotonic eigenval- 
ues, we report the result of the calculation at largest q, 
usually q = 5. For instance, for the 0^~ state at unitar- 
ity, the q = 5 value is E = 7.036 hu, and appears to be 
an upper bound based on the direction of convergence. 
This compares favorably with the result E = 7.010 hw 
(a 0.4% difference) of Ref. 0, in which the center of 
mass coordinate is separated out and a stochastic vari- 
ational approach is employed. Note we add a center of 
mass excitation of 1.5 fio; to the results of Ref. [l^ to allow 
comparison. 

An example of a fast, slightly non-monotonic energy is 
the 1 q state, which was calculated in Ref. to have a 
value E = 6.588(20) hot. Our q = 5 estimate of the new 
interaction is 6.582 frw, well within error. 

Compared with the renormalized contact interaction, 
the new interaction provides much improved convergence 
for four particles, as it did for the three particle system. 
For example, Fig.[3]shows the convergence of the ly state 
using the new and the renormalized contact interactions 
on either side of the crossover. As with all the energies 
we have studied, this displays a marked improvement in 
convergence. On the BEC side, in particular, we find a 
7% difference between the two interactions at q = 4. 

In Fig. 2] we present results for a selected set of low- 
lying energies of the four-particle system as a function of 
inverse scattering length. As for the three-particle case, 
we find crossings between levels with different quantum 
numbers, for instance, the if state and the 0^ state at 
approximately a osc /a = 0.36. In contrast to the three- 
particle system, the angular momentum and parity of 
the ground state (here 0q) remain fixed through the full 
range of scattering lengths. 
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FIG. 3. Convergence of the lowest 1 + state for the four par- 
ticle system. Solid circles are the results calculated with the 
effective interaction, and open circles are the results from 
the renormalized contact interaction, (a): Convergence of 
the energy vs q on the BCS side of the crossover at 

«osc /a = —0.8. (b): The BEC side of the crossover at 
aosc/a, = +0.8. Inset: The absolute values of the differences 
|A£ (9) |. In both cases, the convergence becomes slightly non- 
monotonic at the q = 5 point. 




FIG. 4. (Color online) Selected low-lying energies of the four- 
particle system vs. scattering length. Symbols are the esti- 
mates for q — ► oo as described in the text. Lines are linear 
interpolations. 



B. Thermodynamics 

As thermodynamics is one of the central topics of inter- 
est for cold atoms, we test the usefulness of the new sep- 
arable effective interaction for thermodynamical studies 
by studying the convergence of the free energy (at a given 
temperature) versus q for the three-particle system. Our 
calculations were done at temperatures between T = hoj 
and T = 2.0 fko as a function of q and N m n X via a direct 
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FIG. 5. Convergence of the free energy in iV max at T — OAhw 
and q = 4 calculated with the effective interaction (solid 
circles) and renormalized contact interactions (open circles). 
The solid lines are fits to the interpolating function log( A_F) = 
A + Bq + Cq 2 with fit parameters A,B,C, and are used to ex- 
trapolate to the limit iV max —¥ oo. 



FIG. 6. Free energy versus temperature for the three-particle 
system. Solid circles: new effective interaction at q = 4. Solid 
line: exact result. Dashed line: non-interacting system. The 
inset demonstrates the convergence of F vs. q at T = 0.4 hw. 
The solid circles were calculated using the separable effective 
interaction, while the open circles are from the renormalized 
contact interaction. The horizontal line indicates the exact 
result of 3.6537 . . . hw. 



diagonalization of the Hamiltonian through N max = 10. 
We compare our calculations with exact results for the 
three-particle system. 

An example of the convergence in iV max for q = 4 and 
T = QAhw is shown in Fig. [5j where we observe smooth 
behavior for both interactions. For higher temperatures, 
the smooth character of the AF vs. -/V max is preserved, 
while the curvature flattens out to a straight line (on a 
logarithmic scale), allowing for accurate extrapolations 

in N max . 

The convergence as a function of q is demonstrated in 
the inset to Fig. [6] for T = 0.4 hw. The convergence 
using the new effective interaction is significantly faster 
than for the renormalized contact interaction. In gen- 
eral, we found that the convergence in q for the new in- 
teraction was quite uniform and fast for temperatures up 
to approximately T = l.Ohw. In fact, at higher tem- 
peratures, the abundance of interactionless states (see 
below) diminishes the error for even small q values so 
that at T = 2.0 hw, the difference between q = 1 and 
the exact value is only about 0.8%, compared to 2.1% at 
T = 0.1 hjj. By comparison, the difference in the free en- 
ergy between the non-interacting and interacting systems 
is 2% at T = 2.0 hw and 30% at T = 0.1 hw. 

Our thermodynamics studies arc summarized in Figs. [6] 
and where we compare our calculations of the free 
energy and entropy as a function of temperature with 
exact calculations for the interacting and non-interacting 
systems (described below). We find excellent agreement. 
Although this is only a three-body system, the entropy 
curve indicates that our results apply in the entirety of 
the interesting quantum degenerate regime. 

Since the three-body problem has been solved exactly 
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FIG. 7. Entropy versus temperature for the three-particle sys- 
tem, calculated via a numerical derivative of the free energy, 
S — —dF/dT. Symbols and lines are as in Fig. [6] 



at unitarity [2Cj, |2l| , we can compare the thermodynamic 
functions calculated using the new effective interaction 
with exact results. 

We summarize the exact calculation as follows. There 
arc two classes of eigenstates for the trapped Fermi gas, 
interacting states and interactionless states. The in- 
teracting states satisfy the Bethe-Peierls contact condi- 
tion [H 

¥>(ri,r 2 ,r 3 ) = ( — --) A(R ij ,r k ) + 0{r, l} ) 
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17.792 


-11.109 


14.451 


18.459 


-10.874 
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TABLE I. A few values of the exact thermodynamic functions 
of the unitary and non-interacting three-particle systems. E, 
F, and S are the energy, free energy and entropy of the unitary 
system; Ef lee , F{ ICC , and Sfrcc are the same quantities for the 
non-interacting system. 

(as Tij — > 0) via a function A which is non-zero. The 
interactionless states satisfy a similar condition but with 
A identically zero. They are also eigenstates of the non- 
interacting Hamiltonian and are identical for all values 
of the scattering length. 

To calculate the energy levels of the interacting states 
we follow the procedure of Ref. HtJ However, as noted 
there, at large temperatures the number of interaction- 
less states grows much larger than that of the interacting 
states, so that energies for both are necessary to calcu- 
late thermodynamic quantities for the system. The wave 
functions and energies of the interactionless states for the 
three-particle system have been characterized completely 
in Ref. [2lT Below we summarize the main results. 

The energy of any state of the system (interacting or 
interactionless) is given by 

E = (7 + n c + 2v + 9/2)Huj 

where n c is the number of oscillator quanta in the 
center-of-mass excitation, 7 is a scaling exponent, and 
v = 0, 1,2, ... is the hypcrradial quantum number. The 
corresponding intrinsic wave functions are of the form 
tp cx 4>L i l 7 +2) (R 2 /a 2 osc )e < -- R2 / 2a2 ^\ where R is the hy- 
perradius. The function <f) (cx R 1 ) satisfies the Bethc- 
Peierls contact conditions as well as additional equations. 
States with total internal angular momentum I = 0, 1, . . . 
have exponents 7i jn where n = 0, 1, . . .. For the interact- 
ing states 7i.„ = si_ n — 7/2 where s; jn are an infinite 
sequence of solutions to an Z-dependent transcendental 
equation obtained from the three-particle Schrodinger 
equation poj ). For the interactionless states, 7z jTl are in- 
tegers > 2 and have degeneracies; these degeneracies are 
non-trivial and must be determined in order to calculate 
thermodynamic quantities for the system. We followed 
the method of Ref. Hi] to determine these degeneracies 
(see, e.g., Table II of Ref. HH). 

Using the above classification of interacting and in- 
teractionless states, we computed the entire three-body 
spectrum to the extent needed to ensure convergence of 
the free energy at any temperature of interest (see Figs. [5] 
and [7]) . Table U provides a few values of the free energy 
F and entropy S at various temperatures for reference. 
For the non-interacting calculations we use simple closed 
formulas that are derived in the Appendix. 



V. CONCLUSION 

We have extended a new separable effective interac- 
tion for the unitary trapped gas in the configuration- 
interaction framework to regions away from unitarity, 
i.e., the BEC-BCS crossover. The main advantage of 
this interaction, when compared with the conventional 
rcnormalized contact interaction, is its much improved 
convergence of the many-particle energies in the regu- 
larization parameter. This allows accurate calculations 
in the configuration interaction approach. In particular, 
we calculated the low-lying energies for three and four 
particles as a function of the inverse scattering length. 

We also studied the three-particle system at finite tem- 
perature in the unitary limit and found similar conver- 
gence properties of the free energy as a function of the 
regularization parameter. This separable effective inter- 
action may therefore facilitate accurate studies of the 
thermodynamics of the trapped cold atom system. 
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Appendix A: Non-interacting thermodynamics 

In Figs. [S] and [7] we compared our calculations with 
the non-interacting three-body trapped systems (dashed 
lines). For this purpose we derived simple formulas 
for the canonical partition functions of small systems 
of trapped non-interacting spin-1/2 particles. For three 
particles (ttl)i we proceed as follows. Label the state 
of each particle as s = (n x ,n y ^n z ,a)^ where n x ,n y ,n z 
are the 3-D harmonic oscillator quantum numbers and 
a" = ±1/2 is the z-component of the spin. Then the 
canonical partition function is 

«1,«2,«3 

<r 1 +a 2 +(T 3 = l/2 
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where E s = (n x + n y + n z + 3/2)hu> and S(si, s 2 , S3) is 
the selection function 



S(s 1 , s 2 , s 3 ) = 



if any two s,; are equal 

1 otherwise 



= II( 1 -^)- 



(Al) 
(A2) 



<<j 



To perform the sum, one expands the selection function 
to obtain eight terms involving sums over the Maxwell- 
Boltzmann factor e ^l 3 ( E "i+ E "2+ E '3) times a product (say 
D) of zero to three Kronecker 5's. By collapsing the sum 
over the <5's one obtains a geometric series which then 
easily yields a closed-form expression. As an example, for 



D = <5si,s 2 > an d with the notation |n.;| = n 



the corresponding sum is 

S1.s2.s3 

0"l+CT2+CT 3 = l/2 

1 



1 - e- 2 ? 



1 - e~0 



Note that in this case, the sum over er 2 arL d 03 is ex- 
actly 1 (and is zero for the two and three-<5 terms for this 
system). 

The final result is 



1 



1-e-P 



1 - e-^ 
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